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MATHEMATICAL  ASPECTS  OF  FINITE  ELEMENT 


METHODS  FOR  INCOMPRESSIBLE  VISCOUS  FLOWS 


Max  D.  Gunzburger 
Carnegie-Mellon  University 
Pittsburgh,  PA  15213 

We  survey  some  mathematical  aspects  of  finite  element  methods  for 
incompressible  viscous  flows,  concentrating  on  the  steady  primitive  variable 
formulation.  We  address  the  discretization  of  a  weak  formulation  of  the  Navier- 
Stokes  equations;  we  then  consider  the  div-stability  condition,  whose 
satisfaction  insures  the  stability  of  the  approximation.  Specific  choices  of 
finite  element  spaces  for  the  velocity  and  pressure  are  then  discussed. 
Finally,  the  connection  between  different  weak  formulations  and  a  variety  of 
boundary  conditions  is  explored. 


We  wish  to  acknowledge  the  support,  both  for  the  preparation  of  this  paper 
and  for  our  work  on  finite  element  methods  for  the  Navier-Stokes  equations,  to 
the  Air  Force  Office  of  Scientific  Research  under  grant  AFOSR-81-0101  and  to 
the  Institute  for  Computer  Applications  in  Science  and  Engineering  under 
contracts  NAS1-17070  and  NAS1-18107.  We  also  wish  to  acknowledge  the  many 
fruitful  interactions  with  colleagues  and  friends,  most  notably  J.  Boland,  M. 
Cayco,  G.  Fix,  R.  Nicolaides  and  J.  Peterson. 


Accession  For 

NTIS  GRA&I 
DTIC  TAB 
Unannounced 
Justification. 


By - 

Dlstribut i on/ 

Availability  Codes 
Avail  and/or 
Dist  Special 


MATHEMATICAL  ASPECTS  OF  FINITE  ELEMENT 


METHODS  FOR  INCOMPRESSIBLE  VISCOUS  FLOWS 

One  of  the  most  successful  and  well  developed  mathematical  theories 
concerning  finite  element  methods  is  that  connected  with  incompressible  flow 
problems.  The  success  of  this  theory  lies  not  only  in  the  accumulated  elegant 
mathematical  results,  but  also  in  its  impact  on  practical  computations.  The 
outstanding  monographs  by  Girault  and  Raviart  [GR1.GR2]  give  a  rigorous  account 
of  this  theory  and  to  this  day  remain  the  definitive  sources. 

In  this  survey  we  examine  certain  mathematical  aspects  of  finite  element 
methods  for  the  approximate  solution  of  incompressible  flow  problems.  Our 
principal  goal  is  to  present  some  of  the  important  mathematical  results  which 
are  relevant  to  practical  computations.  In  so  doing  we  also  discuss  useful 
algorithms.  Due  to  space  limitations  we  focus  on  the  steady  primitive  variable 
formulation.  Moreover,  even  within  this  narrow  context,  we  will  concentrate  on 
only  one  of  the  very  many  different  known  approaches.  Some  other  approaches  are 
discussed  in,  e.g.,  [GR1 ,GR2,To] . 

We  state  at  the  outset  that  we  make  no  attempt  at  being  comprehensive  in 
our  coverage  or  in  our  attributions.  To  anyone  who  takes  offense,  we  sincerely 
apologize . 

I  -  The  Primitive  Variable  Formulation 

Let  G  denote  a  bounded,  possibly  multiply  connected,  domain  in  R*.  <i-2  or 


3,  and  let  T  denote  Its  boundary.  As  a  prototype  for  incompressible  flow 
problems  we  consider  the  Navier-Stokes  equations 


(1.1)  u  •  grad  u  *  grad  p  =  v  Au  *  f  in  (3 

together  with  the  incompressibility  constraint 

(1.2>  div  u  =  0  in  (i 

and  the  boundary  condition 

<1.3>  u  =  0  on  T 

where  u  is  the  velocity  field,  p  the  pressure,  f  the  given  body  force,  and  u 
the  given  constant  kinematic  viscosity.  In  (1.1)  the  constant  density  has  been 
absorbed  into  the  pressure.  Whenever  u  and  p  represent  non-dimensionalized 
variables,  then  v  is  the  inverse  of  the  Reynolds  number  Re. 

Following  our  detailed  discussion  of  the  approximation  of  solutions  of 

(1.1) — (1.3)  by  finite  element  methods,  we  will  consider  other  incompressible 
flow  formulations,  especially  as  they  concern  boundary  conditions  other  than 
(1.3). 

1.1  ~  Janet  Lon  spacer,  norms  and  (.arms 

In  order  to  introduce  a  Galerkin  type  weak  formulation  through  which  a 
finite  element  approximation  is  determined,  we  first  need  to  define  some 
function  spaces,  associated  norms  and  forms  involving  functions  belonging  to 
those  spaces.  Lucid  and  more  detailed  accounts  concerning  these  spaces  may  be 


-4- 


found  in,  e.g.,  CBaC,GR2,0Rl . 

2 

First  we  denote  by  L  ffi)  the  space  of  functions  which  are  square  integrable 
over  fl  and  which  is  equipped  with  the  inner  product  and  norm 

<p,q)  =  J  pq  and  •  qi! 0  =  (q,q)1/Z  , 

Ci 

respectively.  We  then  define  the  constrained  space 

L2n(C2)  =  [  q  e  LZ(f3)  I  f  q  =  0  )  . 

01  Q 

2 

Thus  Lg(fl)  consists  of  square  integrable  functions  with  zero  mean  over  Q.  This 
space  is  used  in  connection  with  the  pressure;  such  a  constraint  is  needed 
since  it  is  clear  from  (i.l)-<1.3>  that  the  pressure  can  be  determined  only  up 
to  an  arbitrary  constant.  Other  constraints,  e.g.,  fixing  the  pressure  at  a 
given  point,  may  be  used  instead  without  effecting  any  appreciable  change  in 
the  results  discussed  below.  Next  we  define  the  Sobolev  spaces 


Hk((2)  =  (  q  e  LZ(Q)  I  DSq  <=  L2<Q)  for  s=l . k  ) 


where  DS  denotes  any  and  all  derivatives  of  order  s.  Thus  h(Q)  consists  of 
square  integrable  functions  all  of  whose  derivatives  of  order  up  to  k  are  also 
square  integrable.  H  <fl)  comes  equiped  with  the  norm 


where  the  summation  extends  over  all  possible  derivatives  of  order  k  or  less. 

0  2  1 
Clearly  H  (G)=L  <Q).  Of  particular  interest  is  t tie  spare  H  mnsisting  of 
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functions  with  one  square  integrable  derivative  and  the  subspace 


=  (  q  e  H*  (Q)  I  q  =  0  on  T  ) 


whose  elements  have  one  square  integrable  derivative  over  ft  and  which  vanish  on 
the  boundary  f.  These  spaces  have  the  associated  norm 


(1.4) 


,3q  2  ^1/2 


i=l  1 


We  note  that  for  functions  belonging  to  H^(ft)  the  semi-norm 


(1.5) 


,4-  3q  2  a  1/2 

'■"i  *  (  E  :ST  o  ) 

i=l  1 


is  actually  a  norm  equivalent  to  (1.4)  and  thus,  for  such  functions,  (1.5)  may 
be  used  Instead  of  <1.4). 

For  vector  valued  functions  we  use  the  spaces 


flNft)  =  CH^ fft> ld  =  (  v  I  v  e  HK(ft)  for  i  =  l . d  ) 


and 


H^<ft)  =  [Hg(ft)]d  =  (  v  I  v  e  Hg(ft)  for  1=1,..., d  )  . 


For  example,  H^(ft)  consists  of  vector  valued  functions  each  of  whose  components 
belongs  to  H^(ft).  H^(ft)  is  equiped  with  the  norm 


■  (  E  Vk ) 


2  ,1/2 


i  =  l 


alternately,  H^(ft>  has  the  norm 
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,  A  2  \1/2 

iv'i  ■  (E  "A ) 


i=l 


?  /)  *7  /4 

Also,  the  Inner  product  for  functions  belonging  to  L  (C3)~H  (O)-Il  (Q)v  Is 


given  by 


(v,v>  *  J  u»v 
0 


where  there  is  no  ambiguity  possible  resulting  from  using  the  same  notation 
both  the  inner  product  of  scalar  and  vector  valued  functions. 

We  now  define  the  bilinear  forms 


for 


(1.6) 

and 

(1.7) 


a(u,v)  =  v  /  gradu:gradv  for  all  u,  v€H  (  SI  ) 
Si  0 


b(v,q)  =  -f  qdivv  for  all  veH^(G)  and  qeL^(P), 
0  00 


and  the  trilinear  form 


(1.8) 


c(w,u,v)  *  J  v-gradvv  for  all  u,  v,  veH^  (3) . 

a  0 


In  (1.6)  and  (1.8)  we  have  that  (gradu)J  =9u  /3x,  and 

ij  j  i 


d  3u  3v  d 

r~  j  J  r —  OU . 

gradu:gradv  =  )  - —  - —  and  v-gradu*v  =  )  w  — 1  v 

•—  3x  3x  l—  jdx  i 

i , j=l  J  J  1,J=1  j 


■/vV 

vW 

?:>>>: 


v\*"vsS 


v'_ 
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f; 
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iS 


y. 


iV«V»' 


Using  the  bilinear  form  b(*,*>,  we  can  define  the  subspace 


Z  =  (  v  e  ffj(fl)  I  b(v,q)  =  0  for  all  q  <r  L*(Q>  ) 

which  consists  of  (weakly)  diwrq&nce,  functions,  i.e.,  functions  whose 

2 

divergence  is  orthogonal  to  all  Lg(fl)  functions.  Certainly  any  divergence  free 
function,  in  the  strong  sense,  belongs  to  Z. 

1.2-4  Salerkin  tppe  veak  formulation 

The  most  commonly  used  weak  formulation  of  (1.1)-(1.3>  is  the  following. 

2  A  2 

Given  feL  ((2),  we  seek  ueHq((2)  and  peLg(fl)  such  that 

(1.9)  a(u,v)  +  c(u,u,v)  +  b(v,p)  =  (f,v>  for  all  veH^tQ) 

(1.10)  b(u,q>  =  0  for  all  qeL^(C2). 

By  virtue  of  (1.10)  we  see  that  the  solution  u  belongs  to  Z. 

2 

We  note  that  L  ((2)  is  not  the  largest  function  space  for  the  data  f  such 

that  the  problem  (1.9)-(1.10)  makes  sense;  indeed,  all  that  is  required  of  the 

data  is  that  the  right  hand  side  of  (1.9)  be  bounded  and  this  is  possible  for 

some  functions  which  are  not  square  integrable.  However,  for  our  purposes, 
2 

feL  <(2)  is  sufficiently  general. 

It  can  be  easily  verified  that  whenever  a  pair  u,p  satisfies  (1.9)-(1.10> 
and  is  sufficiently  smooth  to  allow  for  the  appropriate  integrations  by  parts, 
then  u,p  is  also  a  solution  of  <1.1)-<1.3>.  Of  course,  (1.9)  —  (1.10)  admit 
solutions  which  are  not  sufficiently  smooth  to  be  solutions  of  (1.1)— (1.3); 
hence  the  terminology  u/eak  formulation  and  general  (fed  solution  are  applied  to 


(1.9)-(1.10)  and  their  solution,  respectively.  On  the  other  hand,  it  is  also 
clear  that  any  solution  of  (i.l)-(1.3),  i.e.,  a  strong.  solution,  satisfies 

<1.9)-<i.i0) . 

For  the  weak  formulation  (i.9)-(1.10) ,  the  boundary  condition  (1.3)  is  an 
essential  one,  i.e.,  it  must  be  imposed  on  the  candidate  solution  functions. 
Below,  in  section  IV. 3,  we  will  discuss  the  natural  boundary  conditions 
associated  with  the  weak  formulation  (i.9)-(i.i0) . 

We  will  not  enter  into  details  concerning  the  existence,  uniqueness, 
continuous  dependence  on  data  and  regularity  of  solutions  of  ( 1 . 9 )-< 1 . 10 ) .  Such 
results  may  be  found  in,  e.g.,  the  definitive  treatise  of  Teman  CTel. 
Furthermore,  many  of  these  results  are  similar  to  those  discussed  below  for  the 
approximate  problem. 


II  -  The  Finite  Element  Problem  and  the  Div^stability  Condition 


II. 1  -  7 he  discrete  finite  element  problem 

Once  the  Galerkin  formulation  (1.9)-(1.10)  is  established,  the  approximate 
problem  which  determines  the  finite  element  solution  is  defined  in  the  usual 
manner.  First  one  chooses  the  approximating  finite  element  spaces,  or  more 
precisely,  a  family  of  finite  element  spaces,  V^1  and  S^1  for  the  velocity  and 
pressure,  respectively.  Here  h  is  a  parameter  which  is  usually  related  to  the 
size  of  the  grid  associated  with  the  finite  element  partitioning  of  0.  Then  one 
requires  that  (1.9)-  (1.10)  hold  for  functions  belonging  to  these  finite 

dimensional  spaces,  i.e.,  one  seeks  u^eV^1  and  p^eS^  such  that 


(2.1) 


h  h 

a (u  , v  ) 


.  .  h  h  h 

♦  c( u  , u  , v  )  + 


.  ,  h  h 
b(v  ,p  ) 


=  (f 


vh> 


for  all 


h 

v  e 


V* 


tv. , 


I 


v  '”5' 


v. 


L™  * 

g 


*3 


Sv 

.>\v 

i:--- 


■W 


<2. 2) 


Mu  ,q  )  =  0 


for  all  q  eS  . 


If  V*1  and  S*1  are  subspaces  of  the  underlying  infinite  dimensional  spaces  of 


(1.9)-  (1.10),  i.e.,  if  V^cH^(Q)  and  S^cLg(Q),  then  the  finite  element  solution 


defined  by  (2.1)-(2.2>  is  said  to  be  conforming.  Otherwise,  i.e.,  if  vVh!  (Q) 


ti  2 

and/or  S  </L^((2),  then  the  method  is  said  to  be  non-conforming.  We  will  restrict 


our  attention  to  examples  of  the  former. 


Once  one  chooses  specific  bases  for  V  and  S  ,  <2.1>-(2.2)  are  equivalent 


to  a  nonlinear  system  of  algebraic  equations.  Indeed,  if  {q  (x>),  J=1,...,J  and 


(v^(x)},  k=l,...,K,  denote  bases  sets  for  Sh  and  V*1,  respectively,  we  may  then 


write 


J  .  K 

p  =  £  «j<yx)  and  u  =Y  Bkvk<x) 


for  some  constants  a.,  j  =  l  ,...,J,  and  8,,  k=l  ,...,K.  Substituting  into  (2.D- 

J  k 


(2.2)  then  yields 


K_  K 

)  a(v,  ,  v. )  8.  +  Y~  c  ( v  ,v,,v,  )  8,  8 
l—  k  l  k  t—  m  k  t  k  n 

k=l  k,m=l 


(2.3) 


♦  Y.  b(vt,qj>  ~  (f'vt)  for 
j=I 
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.>1  _t:L 


<2. 4) 


Y_  b(vk,qi>  BR  =  0  for  1=1,..., J. 
k=l 


which  constitute  a  nonlinear  algebraic,  in  fact,  quadratic,  system  of  J+K 


equations  for  the  J+K  unknowns  a.,  j=l,...,J,  and  8  ,  k=l,...,K.  Note  that  the 

J  « 


the  discrete  continuity  equation  (2.2)  yields  the  JxK  rectangular  linear  system 
(2.4) . 


II. 2  -  The  ditr-  stab  it  ittf  cand.it  ion 

In  the  positive  definite  case,  e.g.,  for  the  equations  of  linear 
elasticity,  the  mere  inclusion  of  the  finite  element  spaces  within  the 
underlying  function  spaces  is  essentially  sufficient  to  assure  that  the 
approximations  are  well  defined  and  are  as  accurate  as  possible  for  the  type  of 


finite  elements  functions  being  used.  Here  the  inclusions  V^cH^(Q>  and  S^cL^(Q) 


are  not  by  themselves  sufficient  to  produce  stable,  meaningful  approximations. 
We  find  ourselves  in  the  realm  of  what  are  known  as  mixed  finite  element 
me  thodA. 

There  are  number  of  conditions  which  the  elements  belonging  to  the  finite 
element  spaces  should  satisfy.  Most  of  them,  e.g.,  the  boundedfiess  of  the 
various  bilinear  and  trilinear  forms,  are  easily  satisfied  by  conforming  finite 
element  spaces.  The  one  condition  which  presents  a  problem  has  the  following 
mathematical  realization: 


h  h 

given  any  q  eS  , 


(2.5) 


sup 

0*vheV 


,  ,  h  h 
b(v  ,q  ) 


r  D ( v  ,q  )  ^ 

.  .  h..  j 


,  h 
2  Y  q  I 


m 


of  linear 

within  the 

r»  ,i: 
*  •  •/ 

that  the 

the  type  of 

h:-; 

m 


i 


K 


•  V  / 

.V.V 

N'.VC' 


\r* 

V* 


where  the  constant  y>0  may  be  chosen  independent  of  h  and  of  the 
particular  choice  of  q^eS*1. 


This  condition  may  be  equivalently  expressed  in  the  form: 


given  any  q^eS^  there  exists  a  non-zero  v^sV^1  such  that 


(2.6)  b(vh,qh)  z  y!!qh!| 


where  the  constant  y >0  may  be  chosen  independent  of  h  and  of  the 
particular  choice  of  q^eS^1. 


Of  course,  for  each  qh  a  different  vh  may  be  used  in  order  to  satisfy  (2.6). 

The  condition  (2.S),  or  equivalently  (2.6),  is  variously  known  as  the 


f.adpf  he  nAkapa-'Ba.lM.  Aka-Use  i  or  the  £3BB  or  the  in(.-Aup  condition,  the  latter 
designation  following  from  the  third  equivalent  form: 


there  exists  a  y>0,  independent  of  h,  such  that 


(2.7) 


.  ,  h  h . 
f  b(v  ,q  )  -, 

inf  sup  — - - —  2  y 

0* qheSh  Otv^V11  11  v  i;  1 ;;  <1  '  0 


We  will  refer  to  any  of  the  equivalent  statements  (2.5)— (2.7)  as  the  condition 
for  liv-Atakil  tty.  Note  that  these  have  nothing  to  do  with  the  non-linearity  of 
the  Naviei — Stokes  equations  and,  in  fact,  the  possible  problems  its 
satisfaction  poses  is  shared  by  the  linear  equations  of  Stokes  flow. 

Associated  with  the  finite  element  spaces  and  S^1  and  the  bilinear  form 


m 


V;-  ••.< 


!>(•,•)  we  have  the  subspace 


=  (  v  e  V^1  I  b(v^,q^)  =  0  for  all  q^  e  }  . 


of  didcrntet.it  dittergjence  (.ree  (.unctiono.  In  general  Z^Vz,  even  when  V^cH^(Q) 
/i  2 

and  S  cL0(n>,  i.e.,  discretely  solenoidal  functions  are  not  necessarly 

solenoidal.  This  is,  of  course,  entirely  analogous  to  the  finite  difference 
case,  e.g.,  a  function  satifying  a  difference  approximation  to  the 
incompressibility  constraint  is  not  in  general  solenoidal.  A  measure  of  the 
"angle”  between  the  spaces  Z*1  and  Z  is  given  by 


(2.8) 


0  =  sup  inf  "z  -  zh:: 
2heZh  zeZ 

II  zhll  =1 
1 


In  general,  O£0*l,  which  is  easily  seen  by  observing  that  for  Z  eZ,  9=0,  and 
that  by  choosing  z=0,  0=1. 

Note  that  because  of  (2.2),  the  approximate  velocity  u^eZ^,  i.e.,  uh  is 
discretely  solenoidal.  However,  since  in  general  zVz,  divu^xO.  Loosly 
speaking,  the  div-stability  condition  (2.5)  ensures,  as  h->0  at  least,  that 
discretely  solenoidal  functions  tend  to  solenoidal  functions. 


II. 3  -  "£rror  edtimated  and  other  reouttd  concern  ing.  the  approximate  oolut  ion 

We  now  present  some  of  the  available  mathematical  results  concerning  the 
solution  u^,ph  of  the  finite  element  problem  (2.1)-<2.2>.  Here  we  addwne  that 
the  chosen  finite  element  spaces  j/1  and  Sh  satisfy  the  div-stability  condition 


(2.5).  Subsequently,  we  will  look  into  the  issue  of  verifying  that  condition. 
The  summary  presented  is  based  on  the  detailed  analysis  found  in  [CR,  JR,  GR1, 


GR2,  and  GP]. 

First  off,  for  any  feL^fQ),  <2.1>-<2.2)  has  a  solution  u^,p^,  provided  that 
the  div-stability  condition  (2.5)  holds.  However,  one  can  prove  that  the 
solution  is  unique  only  for  "sufficiently  small"  data  f  or  "sufficiently  large” 
viscosity  v.  More  precisely,  let 


k  = 


sup 

h  h  h  trh 
u  , v  , w  ev 


.  h  h  h 
a(w  , u  , v  ) 

I  uh  I  .  I  vh  I  .  I  w^1 1  . 
Ill 


For  standard  choices  of  finite  element  spaces  <  can  be  shown  to  be  independent 
of  h  and,  in  fact,  depends  only  on  Qc Rd  and  d.  Then,  one  can  show  that  (2.1)- 
(2.2)  has  a  unique  solution  whenever 


K_ 

2 


v 


sup 

h  .  h 
v  eV 


If 

n 


•  V 


s  1  . 


This  condition  is  very  similar  to  the  one  which  is  needed  to  show  the 
uniqueness  of  the  solution  of  <1.9)-(1.10)  and  in  fact  the  latter  implies  the 
former,  i.e.,  whenever  <1.9)  —  <1.10)  can  be  shown  to  have  a  unique  solution, 
then,  provided  the  div-  stability  condition  is  satisfied,  <2.1)-<2.2>  also  have 
a  unique  solution. 

When  one  can  show  that  <1.9>-<1.10)  has  a  unique  solution,  it  can  also  be 
shown  that  the  finite  element  solution  of  ( 2 . 1 ) — ( 2 . 2 )  converges  to  that 
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solution.  In  addition,  something  can  be  said  about  the  convergence  of  the 
finite  element  solution  even  when  (1.9)-(1.10)  does  not  posses  a  unique 
solution.  For  details,  see  CGR1  and  GR21. 

Error  estimates  can  also  be  derived.  Provided  that  the  div-st ability 
condition  is  satisfied,  we  have  that 

(2.9)  Hu  -  uh!l.  s.  C.  inf  ilu  -  v^li  .  +  C_9  inf  Ip  -  q^ll. 

1  1  h  _h  12  f,  h  o 

v  ev  q  eS 

and 

(2.10)  ;lp  -  *  C_  inf  Ilu  -  vh",  +  C  inf  p  -  qh 

u  J  h  1  q  h  -h  u 

v  ev  q  eS 

where  9  is  defined  in  (2.8)  and  C^,  i=l,..,4,  are  constants  independent  of  h. 

These  estimates  are  optimal  for  the  "graph  norm"  'lull  ^-Hpl  Q  of  functions 

belonging  to  H^(Q)*Lg(n>  in  the  sense  that  the  rate  of  convergence  of  the 

finite  element  solution,  measured  in  this  norm,  is  the  same  as  that  of  the  best 

approximation  to  u  and  p  out  of  and  S^,  respectively. 

If  the  solution  of  ( 1 . 9>-(  1 . 10 ) ,  or  more  precisely,  of  the*  linearized 

adjoint  problem  corresponding  to  (1.3)— (1.10) ,  is  sufficiently  regular,  then 

2 

one  can  obtain  an  improved  velocity  error  estimate  in  the  L  (£l)-norm,  namely 

(2.11)  lu  -  uh!l Q  s  C^hilu  -  u^l 
where  again  C,.  is  independent  of  h. 

We  see  that  once  the  div-stability  condition  is  satisfied,  the  error  in  the 
finite  element  approximation  depends  only  on  the  ability  to  approximate  in  the 


•‘.i 


*\  <r.  «\ 


chosen  finite  element  subspaces.  In  general,  (2.9)-<2.10>  indicate  that  the 

velocity  and  pressure  errors  are  coupled.  Furthermore,  one  finds  that  it  is 

efficient  to  equilibrate  the  rates  of  convergence  of  the  two  terms  on  the  right 
hand  side  of  <2.9)-<2.10> .  For  this  reason,  one  would  like  to  use,  for  example, 
polynomials  of  one  degree  higher  for  the  velocity  components  than  those  used 
for  the  pressure.  As  a  final  comment,  we  note  that  the  constants  appearing  in 
(2.9>-(2.10>  are  in  general  proportional  to  1/y  where  y  is  the  stability 

constant  appearing  in  (2.5). 

II. 4  -  Ver  i(.]f-inq.  the  d  itr~  $tah  ii  it  >f  condition 

For  particular  choices  of  and  S^1,  it  is  usually  not  an  easy  matter  to 
verify  that  the  div-stability  condition  holds.  To  accomplish  this  task  for 
families  of  such  spaces  is  even  more  difficult.  Here,  we  sketch  three 

techniques  for  verifying  the  div-stability  condition. 

a)  Fortin's  method  -  One  seemingly  attractive  method  of  showing  that  the 
div-  stability  condition  holds  is  due  to  Fortin.  He  has  shown  [FI  that  the  div- 
stability  condition  (2.5)  is  equivalent  to  the  existence  of  a  linear  operator 
from  H^(£2)->V^  such  that  given  any  veH^(fl) 

b(fl^v,qh)  =  b(v,qfl)  for  all  q^eS^1 

and 

III  vl^  i  Civ1 1 

where  the  constant  C>0  may  be  chosen  independent  of  h  and  of  the  particular 
choice  of  veH^tQ).  Thus  the  task  of  verifying  the  div-stabi 1 lty  condition  (2.5) 
is  reduced  to  the  task  of  showing  the  existence  of  the  operator  fT^ ; 
unfortunately,  although  the  latter  task  has  been  accomplished  m  a  f>>v  specif'. 


settings,  in  general,  it  is  also  a  difficult  thing  to  do. 


(2. IS) 


for  k=0,i  , 


lirw  I.  s  CAh  li/I, 

k  4  1 


to  yield 


(2.16) 


sup 

0*vhe 


,  ,  h  h 
b(v  ,q  > 

,  h. 

Iv  I, 


2  C5  '  CGhl<*  \ 


h  „h 


h. 


for  all  q  eS  with  q  q-1 


VerfUrth  then  shows  that  the  div-stability  condition  (2.S)  follows  from  (2.14) 
and  (2.16)  provided  the  constants  C,...,C  are  Independent  of  h. 

1  O 

Thus  the  main  task  of  applying  his  method,  once  the  inverse  inequality 
(2.12)  and  the  approximation  theoretic  result  (2. IS)  have  been  shown  to  hold 
for  the  discrete  velocity  space  is  to  show  that  (2.13)  is  valid. 

c)  The  Boland-Nlcolaides  method  -  A  more  useful  method,  in  the  sense  of 
having  wide  applicability  and  relative  ease  of  use,  has  been  developed  by 
Boland  and  Nicolaides  [BN11.  One  difficulty  with  verifying  the  div-stability 
condition  (2.5)  is  its  g.tot>a l  nature;  Boland  and  Nicolaides  have  shown  how  to 
locahfe  the  difficult  part  of  the  verification  process. 

Specifically,  consider  a  subdivision  of  (1  into  disjoint  macro- elements  0 
r=l,...,R,  each  of  which  consists  of  one  or  a  few  elements  in  the*  the  finite 
element  triangulation  associated  with  V^1  and  S^.  The  number  of  elements  within 
a  macro-  element  is  independent  of  h,  i.e.,  as  we  refine  the  mesh  the  macro¬ 
elements  are  also  refined  so  that  they  always  contain  the  same  number  of 
elements.  Let  denote  the  boundary  of  the  macro-element  (1 . 

Now,  first  suppose  that  the  div-stabi  1  i ty  condition  holds  for  the  pair  V*1 
and  locatt?  over  a  macro-element,  i.e.,  there  exists  a  constant  r>0, 

independent  of  h  and  of  the  particular  choice  of  macro-element,  such  that 


v\.v\N 


Vvl 


Z  ji 

V  v  ", 


(2.17) 


.  ,  h  h , 

f  b(v  ,q  )  )  h..  .  h  _h 

sup  I  - - -  I  *  y:iq  >Iq  for  all  q  eS^ 

o*vhe^  :;v  r 


where 


vj  =  (  !  v=0  on  Tr  )  and  =  (  qeSh|Q  ) 


Since  and  have  fixed  small  dimension,  independent  of  h,  <2.17)  may  often 


be  verified  by  a  direct  computation. 


Second,  suppose  that  the  div-stability  condition  holds  ilolxilly  for  the 


spaces  V  and  S  where 


(2.18) 


}  M 
l  sh  -  (  Lo 


(Q>  piecewise  constant  functions  with  respect  1 
to  the  macro-elements  Q  ,  r~l . R  j 


i.e.,  suppose  that  there  exists  a  constant  ?>0,  independent  of  h,  such  that 


(2.19) 


sup 

0/vheVh 


,  .  h  h 
b(v  , q  > 


~  h  c  . h  ~h 
t  7  q  q  for  all  q  fS 


Sumarizing  the  Boland-Nicolaides  method,  suppose  we  know  that  the  pair 


is  tocaltp  with  constant  r  independent  of  h,  i.e.,  in  the 


sense  of  (2.17).  Further,  suppose  that  the  compos tson  spaces  V^.S^,  which 


satisfy  (2.18),  are  flo/rally  div-stable  with  constant  ?  independent  of  h,  i.e., 


in  the  sense  of  (2.19).  Then  the  spaces  V^,S^  are  /to6oZZ(f  i  with  a 


ronstant  y  Independent  of  h.  Thus,  through  the  use  of  the  comparison  spaces  the 


div-stability  of  the  pair  V^,S^  need  only  be  checked  locally,  i.e.,  over  a 


macro-element . 


.  •  .t.  .  •  .  - 


t  V* 
v» 

•  V 

V 


'V 

r  'v 

v 

fXt 


s's' 

:'v^ 


■V  \  <• 

v. 


A/C-c. 

A.V 

vc>>d 

4 


*  ~  *  •  *  * 


>>> 


V> 
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Thls  method  has  been  succesfully  used,  e.g.,  in  [BC,BN1 ,BN3,GR21 ,  to  show 
the  div-stability  of  a  variety  of  well  known  elements  and  some  novel  ones  as 
well,  both  in  two  and  three  dimensions. 

II. S  -  Example  a  of  unstable  spaces  including,  the  L>  it  inea/'-t  or.st  ar.t  pair 

There  are  different  ways  in  which  arbitrarily  chosen  finite  element  spaces 

may  fail  to  satisfy  the  div-stability  condition.  Here  we  discuss  some  of  these 

and  then  give  specific  examples,  focusing  on  the  much  studied  and  much 

misunderstood  bilinear  velocity-constant  pressure  pair. 

The  most  catastrophic  type  of  failure  is  for  <2. 2),  or  equivalent lv  <2.4>. 

to  imply  that  u^=0,  i.e.,  the  only  discretely  solenoidal  field  belonging  to  V1 

is  the  zero  vector.  The  approximate  solution  is  useless  since,  of  ''nurse,  u^1  0 

cannot  be  a  good  approximate  solution  of  the  Navler-Stokes  equations.  This  type 

of  situation  can  usually  be  detected  by  a  counting  argument,  i.e.,  the  discrete 

divergence  matrix  b(v.,q  >,  j=l,...,J  and  k=l . K,  appearing  in  <2.4>  has 

"  J 

more  independent  rows  than  columns. 

Less  catastrophic  is  the  situation  wherein  fur  one  or  a  feu,  hut  no!  » I  . 
q^eS^  we  have  that  b<v^,q ^>-0  for  all  so  that  y  0  in  <  2 .  S  >  .  This  kind  of 

failure  of  the  dlv-stabillty  condition  is  usually  easy  to  detect  since  ir 
results,  in  practice,  in  the  discrete  divergence  matrix  being  rank  deficient 
Furthermore,  if  these  type  of  pressure  modes  q^1  are  the  sole  reason  f  ,r  ' 
invalidity  of  (2.5),  one  may  often  still  obtain,  through  a  filtering  pro.es, 
useful  approximations. 

A  more  subtle  failure  of  the  div-stability  condition  is  the  -  ase  vh<-re  f  .r 


at  least  some  q  rS 


(2.20) 


,  .  h  h . 

h  r  b(v  ,q  >  •.  h 

C  h'iq  1  *  sup  h—  )  ‘  CzMq  1Q 

o/ver  iiv  n1 


for  some  constants  and  independent  of  h.  In  this  case  r-0<h>  where  r  is 


the  constant  appearing  in  (2.5).  In  practice  this  may  result  in  a  loss  of 


accuracy,  especially  for  the  pressure  approximations.  Such  instabilities  are 


harder  to  detect  because,  of  course,  computations  are  usually  carried  out  using 


a  finite  value  of  h.  In  particular  no  problems  such  as  those  caused  by  rank 


deficient  approximations  to  the  continuity  equation  are  encountered.  This  type 


of  situation  points  out  the  dangers  of  calculating  on  only  one  grid  and  of  not 


at  least  performing  serious  mesh  refinement  studies.  It  also  points  out  the 


usefulness  of  rigorous  results  concerning  the  stability,  or  lack  thereof,  of 


finite  element  spaces. 


a)  -  An  unstable  linear-constant  pair  -  An  example  of  the  first  and  most 


catastrophic  instability  is  the  following  seemingly  natural  choice  for  the 


velocity  and  pressure  finite  element  spaces.  Let  0  be  a  square  which  is 


triangulated  as  in  the  figure  below.  For  the  velocity  approximations  we 


piecewise  linear  functions  with  respect  to  the 


given  t r i angulat ion  which  are  continuous  over  Q 


and  which  vanish  on  f.  For  the  discrete  pressures 


we  choose  piecewise  constant  functions  with 


respect  to  the  same  t r l angulat ion  and  having  zer< 


in  over  Q.  Tlearlv  f  1}  )  and  S^<  .  For 


tills  choice  the  onlv  discrete  veloiltv  field 


h  h 

u  •  v  satisfying  the  discrete  ;  n<  o 


rnpt  ess  l  III  1  It  V 


i  fist  r  a  l  rif  '2.2'  is  u  0 , 


f  fit-  Hi  \ 


r  ef  e  1  V  so  1  e(|i  I  1  d  i  .  t  V  f  1  1  1  .  S  *  he  .’•■! 


_ ^  -  ^ 


vector!  One  easily  sees  that  if  there  are  N  cells  to  side,  that  the  number  of 
equations  in  <2. 4)  is  J=dim(S^>  =2N^  -1  which  is  greater  than  the  number  of 
columns  K=dim( V*1  )=2(N-1  )^. 

In  the  above  example  we  see  that  the  discrete  incompressibilty  condition 
(2.2)  imposes  too  many  constraints  relative  to  the  available  velocity  degrees 
of  freedom.  In  fact,  dim(S^ > >dim( ) .  In  order  to  remedy  the  situation  one 
must,  at  least,  increase  the  dimension  of  V*1  relative  to  that  of  s\ 

b)  The  blllneai — constant  element  pair  -  We  next  consider  the  bilinear 
velocity-  constant  pressure  pair  which  is  often  refered  to  as  the  Q^-Pq  element 
pair.  Again  consider  the  case  of  12  being  a  square  and  consider  the 
"triangulation"  of  the  figure  below.  We  now  choose  to  consist  of  piecewise 
bilinear  functions  with  respect  to  this  triangulation  which  are  continuous  over 
12  and  which  vanish  on  T.  For  we  choose  piecewise  constant  functions  over 
the  same  triangulation  and  which  have  zero  mean  over  (2.  Once  again  the 
inclusions  V^c H^(fl)  and  ShcL^(W)  hold.  The  simple  counting  argument  used  for 

the  first  example  does  not  yield  any  definitive  information  since  dimtv  )= 

2  h  2 

2 ■  N- 1  >  ,  ‘  *.►•  -rune  as  before,  while  now  dim  >'5  >-N  -1. 

It  is  well  known,  e.g.,  see  I F ,  BH,  SGLGE, 

IP,  GNP 1 ,  that  this  bilinear-constant  element  , 

pair  exhibits  the  disastrous  "checkerboard"  mode. 


1  .  M  . 

,  for 

t  he 

part  : 

tcular  disci 

■et  e 

press 

ure  f 

leld 

h  c- 

q  f5 

^  wti  1  • 

'  h  13 

*1  i 

n  the  "red 

squares" 

and  - 

1  i  n 

t  he 

"blac 

k  squares" 

we  have  t 

hat 

b<v^. 

qh)  0 

for 

all 

VheVh 

This 

;  is 

an  example 

of 

t  tie  S  e 

t  olid 

type 

>f  : 

nst  at) 

1  1  1 1  v 

<i  l  sc 

ussed  above 

The  sit 

igle 

bad" 

pr  e-; 

SMI  ** 

mode 

»  an 

be  easily 

f  1  1 

t.  ered 

out  , 

aii-i 

*  Ner 

m  f  * » r  m 

some 

have 

suggest  t-<l 

r  hat 

o|t<  e 

tins 

-22- 


is  taken  care  of,  the  bilinear-constant  element  pair  can  be  safely  used. 

However,  this  is  not  the  whole  story  for  the  bilinear-constant  element 
pair.  Boland  and  Nicolaides  IBN21  have  shown  that  there  exist  other  pressure 
modes  for  which  (2.20)  is  satisfied.  The  left  hand  inequality  of  (2.4)  was 
previously  known  CJP1,  at  least  in  the  different  context  of  penalty  methods.  Of 
course,  the  left  inequality  does  not  imply  the  right,  and  certainly  doesn't 
imply  that  for  those  modes  the  stability  constant  r=0(h>.  However,  Boland  and 
Nicolaides  have  shown  that  this  is  Indeed  the  case.  Moreover,  they  have  shown 
CBN31  that  there  exist  data  f  for  which  the  pressure  approximations  do  not 
converge  and  that  it  is  also  possible  to  set  up  problems  for  which  the  velocity 
approximations  do  not  converge  as  well.  At  the  least,  since  the  constants  in 
the  error  estimates  (2.9>-<2.11>  are  proportional  to  r  *,  there  will  likely  be 
a  loss  of  accuracy  due  to  these  pressure  modes.  Their  conclusions  are  worth 
noting,  especially  in  view  of  the  fact  that  the  bilinear-  constant  element 
pair,  with  the  checkerboard  mode  filtered  out,  has  been  used  on  numerous 
ocassions  in  "practical"  computations. 


Ill  -  Finite  Element  Spaces  for  the  Primitive  Variable  Formulation 

In  this  section  we  discuss  pressure  and  velocity  finite  element  spaces 
which  have  been  rigorously  shown  to  satisfy  the  div-stabil lty  condition.  There 
are  many  such  pairs  known,  especially  for  two  dimensional  problems;  therefore 
we  will  restrict  our  attention  to  pairs  which  have  proven  to  be  of  the  most 
practical  utility. 

Throughout,  P  <T>  denotes  the  space  of  polynomials  of  degree  less  than 
equal  to  k  with  respect  to  the  set  Ir and  IP  (I>1*  denotes  the  space  of  d- 

K 


-23- 


vector  valued  functions  each  of  whose  components  belong  to  P  (3D).  Analogous 

K 


definitions  hold  for  Q^(3D)  and  tQ^(S)]  in  the  case  of  functions  which  are 


polynomials  of  degree  less  than  or  equal  to  k  in  each  of  the  coordinate 


directions,  e.g.,  Q^(D)  denotes  piecewise  bilinear  functions  with  respect  to 


k d 


the  set  3D.  Likewise  we  define  the  spaces  C  (3D)  and  CC  (3D)1  of  k-times 


continuously  differentiable  functions  with  respect  to  the  set  3D. 


For  the  most  part,  the  results  below  hold  for  polygonal  domains  in  R  and 


polyhedral  domains  in  R"  .  Through  the  use  of,  e.g.,  isoparametric  elements, 


they  will  also  hold  for  domains  with  curved  bouhdanes  provided  the  latter 


satisfy  the  usual  smoothness  criteria.  Furthermore,  we  assume  that  all 


subdivisions  of  Q  into  finite  elements  which  are  employed  below  satisfy  the 


standard  .  condit ions .  For  details  concerning  these  issues,  one  may  consult, 


e.g.,  CCi ] . 


III.l  -  Tieceu/ise  linear  and  bilineaj'  velocity  f.  ieldo 


We  begin  with  some  examples  involving  piecewise  linear  or  bilinear  velocity 


fields  with  respect  to  a  subdivision  of  Q  into  triangles  or  rectangles, 


respectively.  In  all  cases  the  discrete  velocity  fields  are  continuous  over  Q. 


In  combination  with  these  type  of  velocity  finite  element  spaces  we  will 


consider  both  discontinuous  piecewise  constant  and  continuous,  over  Q, 


piecewise  linear  pressure  fields.  Every  element  pair  listed  satisfies  the  div- 


stability  condition  (2.5).  Moreover,  provided  the  solution  u,p  of  < 1 . 9>-< 1 . 10 ) 


2  1  12 
satisfies  ueH  (ffJfflfg  <Q)  and  peH  <Q>rtL^<Q),  the  following  error  estimates  for 


the  discrete  solution  u  ,p  of  <2.1)-(2.2>  hold  uniformily  in  h: 


(3.1) 


0(h2  ) 


l  h 

P  ~  P 


.''2' 


Thus,  these  elements  yield  first  order  accurate  pressure  approximations  and 
second  order  accurate  velocity  approximations. 

a)  Piecewise  constant  pressures  I  -  For  the  lineai — constant  element  pair 
mentioned  in  section  II.  5  the  discrete  continuity  equation  overconstrained  the 
approximate  velocity  field.  However,  by  employing  different  grids  for  the 
pressure  and  velocity  fields,  the  linear-constant  element  pair  may  be  made 
stable.  For  example,  consider  a  given  triangulation  7 of  a  polygonal  domain  Q 
into  triangles.  Then  divide  each  triangle  in  7 into  four  triangles  by  joining 
the  midsides,  thus  defining  a  refined  triangulation  7 '  An  example  is 
provided  in  the  figure  below. 


A  pressure  triangle  in  7^ 


Now  define 


The  four  associated  velocity 
triangles  in  f ^ 


Sh  =  (  q  !  q eP0<A>  ,  Ae7h  ;  f  q  =  0  ) 

a 

V*  =  (  v  I  veCPjfd)]2  ,  Ae7h/2  ;  velC°<Q)lZ  ;  v=0  on  f  ] 


so  that  the  pressure  is  sought  among  piecewise  constants  with  respect  to  the 


triangulation  7 ^  while  the  velocity  is  sought  among  continuous  piecewise  linear 
fields  with  respect  to  the  finer  triangulation  7  ^/2'  The  Pair  of  finite  element 
spaces  defined  by  (3.2)  are  known  to  satisfy  the  div-stability  condition  <2.S> 
and  thus  yield  optimally  accurate  approximations  satisfying  (3.1) 

b)  Piecewise  constant  pressures  II  -  For  the  unstable  linear-constant 
element  pair  of  section  II.  S  there  was  one  velocity  element  for  each  pressure 
element;  for  the  stable  lineal — constant  element  pair  (3.2)  there  are  four 
velocity  triangles  for  each  pressure  triangle.  Stable  lineal — constant  element 
pairs  may  be  defined  wherein  the  ratio  of  discrete  pressures  to  velocities  is 
not  so  high.  For  example,  let  the  velocity  space  be  as  in  (3.2);  now  define 
the  pressure  space  S ^  through  the  following  choice  of  basis.  For  each  triangle 
of  7^  we  define  three  basis  functions,  namely  piecewise  constants  which  are 
unity  in  the  shaded  areas  in  figure  below  and  zero  in  the  unshaded  areas.  Of 
course,  outside  the  particular  triangle  of  7^,  the  basis  functions  vanish  as 
well.  This  pressure  space  consists  of  three  out  of  the 


four  possible  piecewise  constants  associated  with  the  four  triangles  in  7 
contained  within  a  single  triangle  in  7^.  Moreover,  there  are  essentially  three 
times  as  many  pressure  of  freedom  for  tins  choice  of 


as  there  are 


for  the  choice  made  in  (3.2).  However,  this  element  pair  is  also  stable,  i.e., 
satisfies  the  div-stability  condition  (2.5)  and  the  error  estimates  (3.D. 

c)  Piecewise  linear  pressures  -  One  may  also  couple  a  piecewise  linear 
velocity  element  with  a  piecewise  linear  pressure  element  and  still  satisfy  the 
div-stability  condition  (2.S)  and  the  estimates  (3.1).  Such  a  pair  was 
introduced  in  IBPj,  analyzed  there  and  in  CV11,  and  is  given  by 

Sh  =  (  q  I  qe?^(A)  ,  ;  qeC°(Q)  ;  J*  q  =  0  ) 

V*1  =  as  in  (3.2). 

Due  to  the  coupling  between  the  pressure  and  velocity  errors  one  cannot  take 

advantage  of  the  better  approximating  ability  of  the  linear  pressure  space. 

Thus,  insofar  as  the  rates  of  convergence,  this  linear-linear  element  pair  is 

no  better  than  the  stable  linear-constant  element  pairs.  However,  in  practical 

calculations  we  have  found  this  to  be  the  best  element  combination  involving 

linear  velocity  fields,  better  in  the  sense  of  giving  more  accuracy  for  useful 

values  of  h.  Furthermore,  this  linear-linear  element  pair  usually  results  in 

fewer  unknowns,  for  the  same  grid,  than  do  the  linear-constant  pairs.  For 

example,  suppose  the  pressure  triangulation  T ^  is  given  by  the  first  figure  of 

section  II. 4  with  N  intervals  on  each  side.  Thus  there  are  2N^  triangles  in  7 

2 

and  the  element  pair  (3.2)  has  2N  -1  pressure  unknowns;  on  the  hand,  the  number 

2 

of  nodes  in  this  triangulation  is  only  (N+l  )  and  thus  the  piecewise  linear 

? 

pressure  space  of  (3.3)  has  only  (N+l)  -1  degrees  of  freedom.  Both  element 

2 

pairs  have  2(2N-1)  velocity  unknowns  so  that  the  linear-linear  element  pair 

2 

(3.3)  has  roughly  N  less  degrees  of  freedom,  for  the  same  grid,  as  does  the 
linear-constant  element  pair  (3.2). 

d )  Piecewise  bilinear  velocity  fields  -  Entirely  analogous  to  rhe 
triangular  elements  described  above,  we  have  the  following  elements  involving 


bilinear  velocity  fields  with  respect  to  rectangular  elements.  More  general 


quadrilateral  elements  may  be  found  from  these  through,  e.g.,  isoparametric 
mappings . 


We  start  with  a  subdividivision  2.  of  Q  into  rectangles,  or  more  generally 

h 


quadrilaterals.  Subsequently  we  divide  each  rectangle  into  four  smaller 
rectangles  by  joining  the  midsides,  thus  creating  another  subdivision  of 

into  rectangles.  See  the  figure  below.  In  all  three  velocity-pressure  element 
pairs 


A  pressure  rectangle  in  2 


The  four  associated  velocity 
rectangles  in  2^/2 


about  to  be  described  we  choose  the  approximating  velocity  space  to  consist  of 
piecwise  bilinear  functions  with  respect  to  the  subdivision  are 
continuous  over  Q  and  which  vanish  on  f,  i.e., 


(3.4) 


=  f  v  1  ve[Q  (0)]2  ,  UeZ.  ;  ve[C°(Q>J2  ;  v=0  on  T 

'  i  n/z 


For  the  first  pressure  space  we  choose  piecewise  constants  with  respect  to 
the  larger  quadrilaterals  of  the  subdivision  2^  and  which  have  zero  mean  over 


E 


WWW 
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Q,  i.e., 


Sh  =  (  q  I  qeQom>  ,  De*h  ;  Jq=0  )  . 

(2 

As  indidicated  in  the  figure  below,  for  the  second  pressure  space  we  choose 
three 


out  of  the  four  possible  piecewise  constants  associated  with  the  rectangles 
belonging  to  2^/2  anc*  ^ave  zero  mean  over  Q.  Finally,  the  third  pressure 

space  consists  of  piecewise  bilinear  functions  with  respect  to  the  subdivision 
2 ^  which  are  continuous  over  G  and  have  zero  mean  over  Q,  i.e., 

(3.5)  Sh  =  (  q  I  qeQ1(0)  ,  0 e2^  ;  qeC°(Q)  ;  Jq=0  ) 

Q 

The  three  velocity-pressure  elements  just  described  satisfy  the  div- 
stability  condition  (2.5)  and  the  error  estimates  (3.1).  Similar  to  the  case 
for  triangles  and  for  the  same  reasons,  the  prefered  element  pair  involving 
bilinear  velocities  is  (3.4)  coupled  with  (3.5),  i.e.,  the  bilinear  velocity- 


bilinear  pressure  pair. 


III. 2  -  7he  Taylor-Hood.  element  paLj « 


We  next  turn  to  quadratic  and  biquadratic  approximate  velocity  fields. 
Suppose  we  have  a  triangulation  7 of  (i.  Then,  the  Taylor-Hood  element  pair 
CTHl  is  defined  by 


<3. 6) 


J  T*1  =  (  v  I  veCP 2fd)]2  ,  AeTh  ;  vg[C(G>]2  ;  v=0  on  f  ) 

1  Sh  =  (  q  I  qeP^(A)  ,  AeX^  ;  q eC((2)  ;  Jq=0  ]  . 

Q 


Note  that  we  are  now  basing  V*1  and  on  the  same  grid  but  on  different  degree 
polynomials,  in  contrast  to  (3.3),  which  uses  the  same  degree  polynomials  but 
different  grids.  The  element  pair  (3.6)  satisfies  the  div-  stability  condition 
<2. 5).  Furthermore,  if  the  solution  (u,p)  of  (1.9>-<1.10>  has  the  indicated 
smoothness,  then  the  following  error  estimates  hold 
uniformily  in  h: 


,.ol  u  -  uh"  j  =  0( h®-1 ) 

1 

,  ueHmrft>nH(i)(ft)  . 

(3.7)  ' 

"u  -  uhl'0  =  0(bm) 

:  whenever  t 

and 

, ,  m-2  or  3 

1 

^  !P  -  Ph:0  =  0(hm)  J 

1 

L  peH®"1  (ft>nL2(ft)  J 

i 

i 

These  results  have  been  established  by  many  authors,  including  [BP.Vl.BNll.  We 
see  from  (3.7)  that  if  ueH2(fl)nH^(fl)  and  peH2  (ft  )DL2  (ft)  then,  in  L2- norms,  we 
have  third  order  accurate  velocity  approximations  and  second  order  accurate 
pressure  approximations.  This  is  an  improvement  over  any  of  the  elements 
involving  linear  velocities. 

One  should  note  that  thetnumber  of  degrees  of  freedom,  both  of  velocity  and 
pressure  type,  associated  with  the  use  of  (3.6)  is  identical  to  that  associated 
with  the  use  of  (3.3),  the  most  efficient  linear  velocity  element.  In  fact,  the 
structure  of  the  discrete  system  resulting  from  a  Taylor-Hood  discretization  is 


in  every  way  identical  to  that  resulting  from  the  use  of  (3.3).  Therefore,  the 


solution  times  for  the  Tayloi — Hood  and  the  linear-linear  discrete  systems  are 
roughly  the  same  if  one  uses  the  same  pressure  triangulation  in  both  cases.  Of 
course,  the  Tayloi — Hood  element  pair  will  yield  better  accuracy  than  the  linear- 
linear  pair,  provided  the  exact  solution  is  sufficiently  smooth. 

On  the  other  hand,  on  the  same  grid,  the  assembly  costs  of  Taylor-Hood  will 
in  general  be  higher  since  one  needs  to  use  higher  order  quadrature  rules  to 
integrate  the  higher  degree  polynomial  integrands  resulting  from  the  Taylor- 
Hood  element  pair.  For  many  solvers,  the  assembly  time  is  overwhelmed  by  the 
solution  time;  therefore  the  increased  assembly  cost  associated  with  < 3 . 6 >  is 
not  a  serious  drawback.  Of  course,  this  is  further  mitigated  by  the  fact  that 
for  the  same  accuracy,  one  may  use  a  coarser  grid  for  (3.6)  than  for  ' 3 . 3  > . 

Summarizing,  provided  the  exact  solution  is  sufficiently  smooth,  the  Taylor- 
Hood  element  pair,  when  compared  to  any  of  the  linear  velocity  elements,  yields 
better  accuracy  for  essentially  the  same  work,  or  alternately,  will  yield  a 
desired  level  of  accuracy  for  less  cost. 

For  rectangles  or  quadrilaterals  we  have  the  analogous  pair 


(3.8) 


/»"■■ 

(  v  1  ve[Qz(D)]2  ,  0s2h 

;  velCin))2  ; 

1  s"  - 

(  q  1  qeQjlD)  ,  Oe.^  ; 

qtC(Q>  ;  j Vo  ) 

Q 

where  2.  denotes  a  subdivision  of  Q  into  rectangles.  This  element  pa.r 
n 

satisfies  the  div-stability  condition  <  2 . 5  >  and  the  error  estimates  <3.?>. 

One  may  well  ask  if  further  efficiencies  may  be  gained  by  using  higher 
order  elements,  e.g.,  cubic  velocities  coupled  with  quadratic  pressures.  Here 
one  needs  to  consider  the  trade-off  between  the  increased  accuracy  of  higher 
order  elements  and  the  increased  complexity  of  those  elements.  \s  in  other 
settings,  e.g.,  structural  mechanics,  one  generally  finds  that  t  fie  optimum 


seems  to  be  achieved  by  quadratic  elements.  Furthermore,  it  is  questionable 
that  in  general  settings  the  exact  solution  of  the  Navier-Stokes  equations  is 
sufficiently  smooth  to  enable  the  potential  better  accuracy  of  higher  order 
elements.  In  our  overall  experience,  we  have  found  the  best  choice  of  velocity- 
pressure  elements  to  be  the  Taylor-Hood  element  pair  or  its 

quadrilateral  counterpart  <3. 8). 


1 1 1. 3  -  D  i  i/«  /tye  nc  e  (.ree  ele^nt^ 

h  h 

Ideally,  one  would  like  to  choose  the  finite  element  spaces  V  and  S  ?<> 
that  the  functions  belonging  to  7^  are  at  least  discretely  divergence  free. 
Certainly  if  the  functions  belonging  to  7^  are  divergence  free  then  they  are 
discretely  divergence  free  as  well,  i.e.,  divv^-0  for  all  implies  that. 

V^-Z^.  Such  a  case  effects  a  great  simplification  since  the  velocity  and 
pressure  uncouple.  Indeed,  we  need  only  solve 


,  h  h  .  h  h  h  .  h ,  f  . ,  h 

afu  ,v  >  +  c< u  ,u  ,v  )  =  (f,v  >  for  all  v  e7 


for  the  discrete  velocit>  ""leld  u^1  since  in  tfus  case  the  term  Mv'Vq1)  in 
<2.1>  vanishes  for  any  v^eV^=Z^.  Also,  since  Z^eZ,  note  that  :  ri  f,he  velocity 
estimate  (2.3),  9=0  so  that  the  velocity  error  depends  only  on  the  ability  to 
approximate  in  v. 

Unfortunately,  although  there  are  known  some  finite  element  pairs  such  that 
the  functions  in  7^  are  at  least  locally  divergence  free,  these  have  proven  to 
be  impractical,  and  we  will  not  consider  them  here.  We  do  mention  that  one 
obvious  method  of  generating  divergence  free  discrete  vector  f  is  t  >  t  ike 

tlie  curl  of  a  piecewise  polynomial  field,  i.e..  ot  i  pie.ewise  poivnomi  1 1 

l ! 


st  reamf  unct  iori .  One  problem  with  this  approach  is 


;  h  it 


■lie 


want  s 


conforming  velocity  field,  i.e.,  >,  then  the  discrete  streamf  unct  ion 

field  must  be  chosen  to  be  conti nuouslv  differentiable  over  Q.  In  J?  this,  of 
course,  necessitates  the  use  of  at  least  quint ic  streamf unct ions  ov*m 
triangles,  or  cubic  polynomials  over  macro-elements,  e.g.,  the  Clough-Toucher 
element.  Non-conforming  velocity  fields  can  also  be  generated  in  this  manner. 


See  [Ca,CN,GRl,  and  GR21  for  details. 


1 1 1  .  A  -  Three  limenjtnnat  e  tew  nt  $ 

Compared  to  the  two  dimensional  setting,  there  are  known  much  fewer  stable 
element  pairs  for  three  dimensional  problems.  However,  there  is  great  interest 
in  this  subject  and  therefore  there  has  been  substantial  recent  progress.  Hero 
we  mention  a  few  of  the  known  stable  three  dimensional  elements. 

In  the  first  place,  the  three  dimensional  analogue  of  tic-  Taylor -Hood 
element  is  known  to  be  stable  in  3-D;  this  may  be  shown  by  the  methods  of 
VerfOrth  or  Boland-Nicolaides .  Specifically,  we  subdivide  »2  into  tetrahedrons 
and  use  continuous  piecewise  quadratic  polynomials  for  the  velocity  and 
continuous  piecewise  linear  polynomials  for  the  pressure.  The  aoor.u  \  . .  f  r}.,.-; 
combination  is  the  same  as  in  the  two  dimensional  case. 

Next  we  consider  linear-constant  elements.  Again,  subdivide  12  inf  > 
tetrahedrons.  For  the  pressure  space  we  choose  piecewise  constants  with 
respect  to  this  initial  subdivision.  Now  we  subdivide  each  tetrahedron  into  ! .? 
smaller  t et r ahedrons  by  first  joining  the  centroid  of  the  faces  t  >  the 
vertices,  and  then  joining  the  centroid  of  the  large  tetrahedron  to  the 
vertices  and  the  centroids  of  t  fie  faces.  For  the  velocity  space  we  choose 
continuous  piecewise  linear  polynomials  with  respect  to  the  smaller 
t et  r ahedrons . 

Another  stable  linear-constant  clement  pa  1  r  is  defined  bv  first  subdividing 


<  4  .  1  > 


7* ( v ,  u,  v  1 


<  W.  V. u  I  J 


“  [  r  ( y ,  U ,  V  )  - 


One  may  easilv  verify  that  <*tu,u,v)  .  <u,u,v>  whenever  tlivu  0  in  12  and  u-  ri  0 
r,  where  n  denotes  the  outward  normal  to  J'  Therefore.  due  t  •  1  . 2  )- <  S  .  3 •  , 
seems  irrelevant  whether  one  uses  «  1 . 8  >  or  •  4 .  1  :  in  i  w«-  iK  t  or  rti...  c 
Navier-  Stokes  equations.  From  an  analysis  point  >t  view.  •  tie  in’ 


<  A .  1  >  is  that  7(ir,u,vi  -f'v.v.u)  for  anv  u.v.weff1 


<12)  while  t  he  Hi  a  logons  :esi)j 


for  (1.8)  holds  on’v  when  divv-0  m  Q  and  one  of  u  0.  v  0  r  w -n  0  •  • .  /' 

We  emphazise  that,  insofar  as  the  a<<  nr  arv  of  t  he  approx  l  mat  :mis 

r  o  nrer  ne  i ,  it  makes  rio  difference  whether  orie  use.-;  t  ♦  8  <  .r  1  .  we  met”! 

point  out  that  mariv  of  the  results  concerning  f  irate  element  approx  l  mat  ions  ■ 
SO  1  lit  ions  of  (  1  .  1  )  -  (  1 . 3  >  we  re  first  oht  a  1  ned  f  hr  oug  ft  t  he  use  if  i  4  1  1  i  in  *  1 

•"  her  hand,  anv  implementation  of  <  4.  1  >  will  result  m  more  .input  ,t  iona.  w«m 
than  the  analogous  implementation  of  <1  8) 


IV.  2  -  f  r.A.nnirt^>  r,e  vu  >  irvlnt  ity  luiur.  ic~ny  <  or.  /  i  f  i  nr,  > 


.  .i*o  ••  it  e  man'. 


•nt  v.  av 


it  .  1 1 f .  i m< 


'otld  1 1  1  otiS  .  Ill  pr  n't  l.  h,  t  he  overwhelming  I  ho  1 1  *■  is  ’o  is e  t  * . •  ■ 
l  rit  e  r  po  1  ant  Ur'  de  sc  r  i  he  this  me  t  hod  for  polvgonal  domains  ,2  K,'  . 


analogous  ideis  ma\  f  >»•  used  •  r.  three  dimensions  and  for  1 . .  m 


i .  • u  ,  t 


rides,  *  he  .  i  ’  t  e  r  t  hi  ' )  i  ig  h  the  aid  of. 


•  par  amet  r  :  •  ••  I  ••  me ; it 


Consider  the  bounder  v  i  oiidd  idi 


u  g  u.  r 


Hid  t  he  set 


V  Uf  H  (..’('  II  ,  at  ,  t 


4  ' 


,1 


Note  that  The  weak  formulation  which  we  will  discretize  is  as 

0  0 

2 

follows:  seek  he  F  and  peL„(Q)  such  hat  <1.9>  and  (1.10)  hold.  Note  that  the 

g  ^  o 

test  function  v  still  belongs  to  ) ,  1 .  e  .  ,  v=0  on  r. 

In  order  to  pose  our  discrete  problem  we  choose  finite  element  spaces 
vW  <Q>  and  S hcL^(fl).  We  denote  by  the  restriction  of  to  the  boundary 

f,  i.e.,  i  consist  of  functions  defined  on  F  arid  which  can  agree  with  the 
boundary  values  of  at  least  one  function  belonging  to  V^1 .  The  finite  element 
functions  belonging  t<>  V^,  being,  for  example,  piecewise  polvnomi  als ,  cannot  m 
general  sat  isf  v  t  lie  boundary  condition  <  4 . 2  >  ;  certainly,  m  general  g,v\. 
Therefore  we  choose  an  approximation  to  g,  which  we  denote  by  g^,  belonging  to 
^ .  The  most  common  choice  for  g^ ,  and  the  one  we  consider  here,  ;s  the 
i  nr  erpo!  ant  of  g  in  l/1 ! 

This  choice  is  trivial  to  implement,  whi-  h  at  least  part i a  1 1 v  account  s  for 

. »  s  ;>opu  l  ar  1 1  v  .  For  example,  suppose  is  a  Lagrange  finite  element  space , 

i.e  one  whose  degrees  of  freedom  are  exclusively  function  values  at  points, 
bet  f  )  ,  k-I,  ,K  denote  tlie  usual  finite  element  basis  for  Let  the  first 

b  .}  these  basis  functions  be  associated  with  inter  ;,.r  nodes  x,  -o  ’fiat  for 

r\ 

FI,  .  k.  v  0  for  »  F  The  remaining  basis  function  Iv  1,  k  K*1 .  .  .  .  K ,  are 

K  K 

issoi  i  at  e<l  with  nodes  1  v  i  rig  on  f  In  practical  implementations  there  are 

more  efficient  node  numbering  schemes  t  hari  the  one  we  are  using,  however  ,  the 
batter  simplifies  ’tie  explanations  being  attempted  here 

t'hoosing  g ^  to  fie  t  fie  tioundar  v  interpolant  of  g  is  t  f,en  e.juivalent  to 

wr  If  ;  rig 


h 

u  x ) 


■  ,  v  'X 

k  k 


>  g  i  x,  1  v,  x 

i_  .  6  k  k 


k  i  ><  k  I  .  .  k  ,  are  r  fie  unknown  <  oef  f  ;  i  i  e  nt  s  *  o  be  dot  .•  r  rn  1  : 


I 


coefficients  of  the  basis  functions  associated  with  boundary  nodes  are  simply 


set  equal  to  g  evaluated  at  the  corresponding  node.  Note  that  (4.4>  implies 
that 


h ,  , 

g  <x) 


K 

Y  *<xk>vk<x)  for  ■ 

k=K+l 


The  contribution  to  u  emanating  from  the  second  summation  of  (4.4)  becomes 

part  of  the  data  of  the  discrete  system  of  equations. 

Once  an  approximation  g^1  is  chosen,  one  may  define  the  set 

V*1  =  (  I  v=g^  on  T  ) 

S 


Note  that  is  the  finite  element  subspace  of  used  in  conjunction  with 

the  homogeneous  boundary  condition  (1.3);  also,  clearly  l^cH1  (Q)  is  not  a 

subset  of  V  .  Now,  the  approximate  problem  may  be  defined  as  follows:  seek 

i/V:v£  and  p^eS^cL^iQ)  such  that  (2.1>~<2.2>  hold  for  all  and  q^eS^, 

respectively.  Again,  the  test  functions  v^1  vanish  on  the  boundary  r. 

The  whole  discussion  of  the  div-stability  condition  (2.5)  carries  over 

intact  to  the  case  of  the  inhomogeneous  boundary  (4.2);  in  ( 2 . 5  >  we  still  use 

the  subspace  of  finite  element  velocity  fields  which  vanish  on  the  boundary. 

Results  analagous  to  those  of  section  IT. 3  can  he  derived  in  a  fairly 

straightforward  manner  with  the  exception  of  some  technicalities  encountered 
2 

for  the  L  (Q)-error  estimate  for  the  velocity  approximation.  Si-e  [CP,  FHP,  GR21 
for  details. 


I n  part irular ,  i f 
then  all  the  results 
spaces  discussed  in 


is  chosen  to  be  the  boundary  interpolant  of  g  in  V^!p, 


e.g.,  error  estimates,  concerning  the  finite  element 
section  III  are  essentially  still  valid  for  the 
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inhomogeneous  velocity  boundary  condition  (4.2).  Again,  see  [GP,  FGP  and  C-R21 
for  details. 

IV.  3  -  Alternate  boundary  cond  it  Lons  and  formutat ionA  of.  the  i/Iacoua  term 

In  this  section  we  examine  how  different  choices  for  the  viscous  term  in 

(1.1)  affect  the  natural  boundary  conditions  of  corresponding  weak  formlations. 
Some  of  this  material  can  be  found  in  CGLN1. 

Due  to  (1.2),  when  v  is  constant,  the  viscous  term  in  <1.1)  may  be  written 
in  the  various  equivalent  forms 

( 4 . 5 . 1 )  vAu  = 

(4.5.2)  divj^  v(  (grad  u)  +  (grad  u>^  ]j  = 

(4.5.3)  -ucurl (curlu)  = 

(4.5.4)  u(  grad(divu)  -  curl(curlu)  )  . 

Although  these  different  realizations  are  equivalent  insofar  as  the  partial 
differential  equations  are  concerned,  we  shall  see  that  each  generates  a 
different  numerical  method. 

If  for  some  reason  v  is  not  constant  or  divu*0,  then  only  (4.5.2)  may  be 
used.  Indeed,  (4.5.2)  is  the  form  of  the  viscous  term  which  arises  naturally 
in  the  derivation  of  the  Navier-Stokes  equations  from  the  principle  of 
conservation  of  linear  momentum  and  the  Cauchy-Poisson  constitutive  equation. 
The  other  three  forms  (4.5.1),  (4.5.3)  and  <4.5.4)  are  derived  from  (4.5.2) 
with  the  aid  of  (1.2)  and  the  assumption  that  inconstant.  In  (1.1)  we  have  used 


tr. 

ft 
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(4.5.1>  only  because  this  is  the  most  popular  choice  in  the  literature;  all  of 
the  results  obtained  so  far  hold  equally  well  if  one  chooses  '4.5.2)  instead. 
As  will  be  seen  from  the  discussion  below,  (4.5.2)  is,  in  general,  to  be 
prefered  to  <4.5.1). 

Denote  two  segments  of  the  boundary  F  by  and  r  .  These  segments  may  be 
empty,  are  not  necessarily  disjoint  and,  in  fact,  may  be  equal.  Now,  for  fixed 
given  functions  g  and  g  ,  define  the  set 

V  =  f  veff1  I  v*n=g  on  F  ;  n*vxn=g  on  f  1 
g  1  °n  n  °x  x  ' 

and  the  spaces 

Vn  =  f  veH1  I  v-n=0  on  T  ;  v*n  =0  on  F  1 
O'-  n  r  i 

and 

S  =  LZn(n>  if  r  =r  ,  S=LZ(fl>  otherwise. 

0  n 

where  v«n  denotes  the  component  of  v  normal  to  the  boundary  r  and 

n*v*n=v-  ( v  n>n  is  the  projection  of  v  onto  the  plane  tangent  to  F.  In  the 

definition  of  V ’  we  may  use  v*n=0  due  to  the  relation  v*n=n*<nw<n),  i.e., 

2 

n*v*n =0  implies  that  n*v=0.  In  R  ,  n*v»n=vr  where  t  is  the  unit  tangent  vector 
to  r. 

Suppose  that  we  wish  to  specify  the  boundary  conditions 

( 4 . 6 . 1  >  u*  n=g  on  r 

°n  n 

and 
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n*u*n=g  on  f 
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i.e.,  the  normal  velocity  on  and  the  tangential  velocity  on  F^, 

respectively.  For  all  the  weak  formulations  which  we  will  consider  involving 
any  of  the  choices  in  (4.5),  (4.6)  will  be  essential  boundary  conditions.  Thus 

the  trial  solution  functions  u  will  satisfy  (4.6),  i.e.,  ue7  ,  and  the  test 

g 

functions  satisfy  yeV^. 

Consider  the  following  weak  formulation:  for  i-1,2,3  or  4,  seek  ueV  and 

S 

peS  such  that 

(4.7)  a.(u,v)  +  b(v,p)  +  c( u,u,v)  =  (f,v)  +  d(v)  for  all  veV^ 
and 

(4.8)  b(u,q>  =  0  for  all  qeS. 


Here,  b(-,-)  and  c(*,-,->  remain  as  in  (1.7)  and  (1.8),  respectively,  and  f 
continues  to  denote  the  body  force  appearing  in  the  momentum  equation.  The 
linear  functional  d(-)  is  given  by 
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(4.10.2)  a^( u,v)  =  1  J*  v (  gradu  +  (gradu)^  ] : [  gradv  +  (gradv)^  ] 

2  n 

(4.10.3)  a^(u,v)  =  uj  (curlu) • (curlv) 

Q 

and 

(4.10.4)  a4(u,v)  =  uj  (curlu)* (curlv)  +  (divu)(divv)  . 

a 

In  the  customary  manner,  should  u  and  p  be  sufficiently  smooth,  one  can, 
through  formal  integration  by  parts  procedures,  ascertain  what  differential 
equation  problem  the  weak  formulation  <4.7)-(4.8)  corresponds  to.  To  begin 
with,  we  know  that  the  boundary  conditions  (4.6)  are  satisfied  since  these  are 
being  required  of  the  candidate  trial  functions  u.  We  also  find  that  the 

differential  equations  (1.1)  and  (1.2)  are  satisfied,  where  in  (1.1)  the 
viscous  term  is  replaced  according  to  (4.5),  depending  on  which  choice  is  made 
in  (4.10).  Finally,  one  finds  the  natural  boundary  conditions  corresponding  to 
the  particular  weak  formulation.  We  will  now  discuss  these  in  some  detail  for 
each  possible  choice  in  (4.10).  • 

Corresponding  to  the  paired  choices  (4.5.1)  and  (4.10.1)  we  have  the 

natural  boundary  conditions 

<4.11.1)  p  -  un*gradu*n  =  r  on  r/r^  and  t>n*gradu*n  -  s  on  T/r^. 

Unfortunately,  these  boundary  conditions  have  no  physical  'waning.  Thus  the 

choice  (4.5.1),  or  equivalently  (4.10.1),  can  only  be  used  in  conjunction  with 

the  boundary  condition  (4.6)  specified  on  all  of  T,  i.e.,  u  given  on  f  =r  =r . 

n  r 


Next,  consider  the  choices  (4.S.1)  and  (4.10.1).  The  natural  boundary- 


conditions  are  then 


(4.11.2) 


|  -p  +  un- (gradu  +  (gradu)^j-n  =  -r  on  r/Tn  and 

1  on- (gradu  +  (gradu)^]*n  =  -s  on  f/r . 


Thus  -r  and  -s  are  the  normal  and  tangential  stresses,  respectively,  on  the 
boundary.  Then,  for  the  choice  (4.10.2),  the  possible  combinations  of  boundary 
conditions  at  a  point  on  the  boundary  f  are  as  follows:  we  may  specify  the 
velocity,  or  we  may  specify  the  normal  velocity  and  the  tangential  stress,  or 
we  may  specify  the  tangential  velocity  and  the  normal  stress.  The  latter 
combinations  are  useful,  e.g.,  for  free  surface  problems  or  at  artificial 
outflow  boundaries.  Details  may  be  found  in  [GLN1 . 

The  third  choice  (4.5.3),  or  <4.10.3),  yields  the  natural  boundary 
conditions 


(4.11.3)  p  =  r  on  r/f  and  ta  =  s/v  on  r/r 

n  r 

so  that  r  and  s  are  the  pressure  p  and  v  times  the  vorticity  a>=curlu, 

* 

respectively,  on  the  boundary.  The  possible  combinations  of  boundary 
conditions  are  now:  we  may  specify  the  velocity,  or  we  may  specify  the  normal 
velocity  and  the  vorticity,  or  we  may  specify  the  tangential  velocity  and  the 
pressure.  The  pressure  is  often  used  as  an  outflow  condition;  the  vorticity  is 
useful  in  exterior  problems  when  matching  to  an  inviscid  irrotational  flow 
since  it  is  well  known  that  the  vorticity  decays  to  its  value  at  infinity 
faster  than  does  the  velocity.  Again,  details  may  be  found  in  fGLNl . 

Unfortunately,  although  the  boundary  conditions  associated  with  the  use  of 

(4.10.3)  can  be  useful,  in  practice  we  cannot  employ  this  particular 


formulation  of  the  viscous  term.  The  reason  for  this  is  that  the  choice 

(4.10.3)  requires  the  use  of  divergence  free  finite  element  velocity  fields  in 
order  for  the  form  to  be  coercive  on  Z*1.  This  condition  is  also  needed 
to  guarantee  the  stability  of  the  approximations  and,  for  the  other  three  cases 
(4.10.1),  (4.10.2)  and  (4.10.4),  is  trivially  satisfied  for  any  choice  of 
conforming  discrete  velocity  space. 

Fortunately,  the  boundary  conditions  (4.11.3)  are  approximately  the  natural 
boundary  conditions  associated  with  the  choice  (4.11.4).  In  fact,  for  (4.10.4), 
we  have  the  natural  boundary  conditions 

(4.11.4)  p  -  vdivu  =  r  on  F/F  and  u  =  s/v  on  r/F 

n  r 

The  second  of  these  is  identical  to  the  second  of  (4.11.3).  If  o  is  "small", 
and/or  if  we  assume  the  incompressibility  constraint  holds  up  to  portions  of 
the  boundary  where  the  normal  velocity  is  not  specified,  then  (p-vdivu)  is 
essentially  equal  to  p.  Thus  we  recover,  at  least  approximately,  the  first 
boundary  condition  of  (4.11.3). 

In  summary,  when  one  has  velocity  and/or  stress  boundary  conditions,  one 
should  use  (4.11.2)  in  (4.7)  and  when  one  has  velocity  and/or  pressure  and/or 
vorticity  boundary  conditions  the  choice  (4.11.4)  is  preferable. 

The  discretization  of  (4.7)-(4.8)  follows  the  usual  procedures  once  one 
chooses  tlie  finite  element  spaces  for  the  velocity  and  the  pressure 
approximations.  The  natural  boundary  conditions  are  automatically  acounted  for 
by  the  inclusion  of  the  linear  functional  d(-)  in  (4.7).  The  essential  boundary 
conditions  on  the  components  of  the  velocity  can  be  enforced  in  a  manner 
analogous  to  that  described  in  section  IV.  2  for  the  case  where  the  complete 
velocity  is  specified  on  the  whole  boundary.  All  material  relating  to  the  div- 


stability  condition  ( 2 . 5 >  is  essentially  still  valid,  and  thus,  insofar  as  that 
condition  is  concerned,  the  particular  choices  of  finite  elements  discussed  in 
section  III  may  still  be  used. 

In  actuality,  there  are  very  few  rigorous  error  estimates  available  for 

boundary  conditions  other  than  the  velocity.  For  polygonal  or  polyhedral 

domains  Q,  the  error  estimates  of  section  II. 3  are  still  valid.  However,  for 

domains  with  c  urved  boundar  ie  a  ,  using  the  type  of  weak  formulations  discussed 

here  may  result  in  a  loss  of  accuracy.  For  example,  for  (4.10.2)  with  normal 

velocity  and  tangential  stress  boundary  conditions,  it  was  shown  by  Verfurth 

[V21  that  there  is  a  loss  of  accuracy  due  to  a  Babuska  type  paradox,  i.e.,  the 

2 

limit  of  solutions  of  problems  posed  on  polygonal  approximations  to  QcR  is  not 
the  solution  of  the  problem  posed  on  0.  Verfurth  IV3I  has  also  shown  how 
through  the  use  of  additional  Lagrange  multipliers  on  the  boundary,  a  different 
weak  formulation  yields  optimal  accuracy. 
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